undef("ispoint_in_polygon")
function is_point_in_polygon(p[2]:float, polygon[*][2]:float, nCount:integer)
local nCross, x, i, p1, p2

begin
    nCross = 0
    x = new(1, "double")
    do i = 0,nCount-1,1
        p1 = polygon(i,:)
        p2 = polygon((i+1)%nCount, :)

        ; calculate the cross point of y=p.y and p1p2
        if (p1(1) .eq. p2(1)) then ; parallel to the X axis
            continue
        end if

        ; the crossover point is on the lower extension line of p1p2
        if (p(1) .lt. min((/p1(1),p2(1)/))) then
            continue
        end if
        ; the crossover point is on the lower extension line of p1p2
        if ((1) .ge. max((/p1(1),p2(1)/))) then
            continue
        end if

        ; calculate the X of crossover point
        x = todouble(p(1) - p1(1))*todouble(p2(0)-p1(0)) / todouble(p2(1)-p1(1)) + p1(0)

        ; only count crossover point in one side
        if (x .gt. p(0)) then
            nCross = nCross + 1
        end if
    end do

    ; if nCross is even number, the point is outside of the polygon
        return ((nCross % 2) .eq. 1)
end

